function [P_S, F_es, wprimeS, educ_mat] = choiceProb_ct_par(Ss, rs, tax, w_mto, mean_wage)
    global a bet0 bet1 xi pe tau w ub lb gamma sigma rts N_w N_a sigmae b phi alpha w_mult pi_e e thetas normalization ss_calc w_bar eta unemp gamma w_gamma vouch_sh
    
    % Pre-allocate matrices
    N_n      = length(Ss);
    u_s      = zeros(N_w, N_a, N_n, N_n, 2);
    educ_mat = zeros(N_w, N_a, N_n, N_n);
    F_es     = zeros(N_w, N_a, N_n, N_n);
    wprimeS  = zeros(N_w, N_a, N_n, N_n);
    ehat     = dot(pi_e, e);
    
    % Sort neighborhoods and calculate variables
    [~, id_max] = sort(-mean_wage);
    [~,id] = min(rs);
    id = id_max(end);
    b0 =  b ;
    w2 = (1 - tax)*w; % Net wage
    wto_share = 0.3;  % Share of income going to MTO housing
    vouch_cap = vouch_sh*rs(id_max(1)); % MTO vocuher cap
    
    % Solve education decision for each potential destination neighborhood

    for c1 = 1:length(Ss)
        for c = 1:length(Ss)
            % Return to unit of education
            E_S = eta * (a' .* (bet0 + bet1*Ss(c)).^xi)  .* f_w(w2) ;         
            % Wage additive constant
            btilde = unemp + b .* f_w(w2) + pe .* E_S;
            
            % Solve for optimal education
            % Quadratic Solution as initial guess
            % MTO: Give the option of moving to the best neighborhood but contribute
            % 30% of income to the rental rate

            if (c1 == id_max(end)) & ((c == id_max(1)))
                
                net_w = w2 - rs(c) + min(max(rs(c) - max( wto_share * w2/(1-tax), rs(c)-vouch_cap), 0), 100*rs(id_max(2))) .* (w <= w_mto);
                discrim = btilde.^2 + (3./tau)*((E_S).^2).*net_w ;
            else
                discrim = btilde.^2 + (3./tau)*((E_S).^2).*(-rs(c) + w2);
            end
            e_init = max(min((-btilde + sqrt(max(discrim,0)))./(3 .* (E_S)) .* (discrim > 0), ub), lb) + 0.1;
          
            % Newtons Method
            for iter = 1:230
               % Precompute
               evar = e_init.^(gamma - 1);
               evar2 = e_init.^gamma;
               %Residual
                if (c1 == id_max(end)) & ((c == id_max(1)))
                    
                    net_w = w2 - rs(c) + min(max(rs(c) - max( wto_share * w2/(1-tax), rs(c)-vouch_cap), 0), 100*rs(id_max(2))) .* (w <= w_mto);
                    foc =  gamma*tau* evar .* btilde + gamma*E_S.*tau .* evar2 - (net_w - tau*evar2) .* E_S * w_gamma; 
                    foc = foc .* (net_w > 0);
                else
                    net_w = w2 - rs(c);
                    foc =  gamma*tau* evar .* btilde + gamma*E_S.*tau .* evar2 - (net_w  - tau*evar2) .* E_S * w_gamma; 
                    foc = foc .* (net_w  > 0);
                end
               foc(isnan(foc)) = 0;

               % Residual derivative
               dfoc =  gamma*(gamma - 1)*tau*e_init.^(gamma - 2) .* btilde + gamma^2*E_S.*tau .* evar  + gamma*tau*evar.* E_S * w_gamma; 
               e_init = (e_init - foc ./ dfoc) .* (net_w  > 0);
               e_init(isnan(e_init )) = 0;

               % Check solution
               if max(abs(foc(:))) < 5e-4
                   break
               elseif iter == 130
                   display(max(abs(foc(:))))
                   error("Did not converge in education")
               end
            end

            e_S = e_init;
            e_S = min(max(e_S, lb), ub);

            % For use in other programs
            F_es(:,:, c1, c)    = (pe + e_S) .* E_S ./  f_w(w2);
            educ_mat(:,:, c1,c) = pe + e_S;
            wprimeS(:,:, c1,c)  = unemp + (b  + F_es(:, :, c1, c)) .* f_w(w2)  ;

            % Calculate utility under each neighborhood choice

            if (c1 == id_max(end)) & ((c == id_max(1)))
                net_w = w2 - rs(c) + min(max(rs(c) - max( wto_share * w2/(1-tax), rs(c)-vouch_cap), 0), 100*rs(id_max(2))) .* (w <= w_mto);
                u_s(:,:, c1, c, 1) = (max((net_w - tau*e_S.^gamma).* ( wprimeS(:,:, c1,c) ), 0*1e-20) );
            else
                u_s(:,:, c1, c, 1) = (max((w2 - rs(c) - tau*e_S.^gamma) .* ( wprimeS(:,:, c1,c) ), 0*1e-20) );
            end
            
            % If wage is too low to afford housing, allocate to cheapest neighborhood
            if (c == length(Ss))
                mysub = sum(u_s(:,:,c1,:,:), 4:5);
                [mrow, mcol] = ind2sub(size(mysub), find(mysub == 0));
                u_s(mrow, mcol, c1, id,1) = 1e-20;
            end
        end
    end
    
    % Amenities shock
    for c = 1:length(Ss)
        u_s(:,:, :,id_max(c), 2) = u_s(:,:, :,id_max(c),1) * thetas(c);
    end

    % Calculate choice probabilities
    my_log_diff = (log(u_s) - max(log(u_s), [], 4)) ; 
    numer = exp(my_log_diff * sigmae  ); %
    P_S = numer ./sum(numer , 4);
      
end